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Abstract 

Ultra-precision machining of metals, the breaking of nanowires under tensile stress and fracture of nanoscale 
materials are examples of technologically important processes which are both extremely difficult and costly 
to investigate experimentally. We describe a multiscale method for the simulation of such systems in which 
the energetically active region is modelled using a robust tight-binding scheme developed at the Naval Research 
Laboratory (NRL-TB) and the rest of the system is treated with molecular dynamics. We present a computer 
code implementing the method, geared towards non-equilibrium, cross-scaled tight-binding and molecular dy- 
namics simulations. Apart from the presentation of the method and implementation, we discuss preliminary 
physical results obtained and discuss their validity. 



1 Introduction 

Ultraprecision machining, or UPM, is a process 
in which a hard (usually diamond or CBN) tool ma- 
chines away a layer of workmaterial as thin as a few 
nanometers. The machined materials include soft 
metals, such as aluminum or copper, glasses and semi- 
conductors. A prime example of a technologically im- 
portant application of UPM is the finishing of alu- 
minum platters for magnetic hard disk drives [1], [2]. 
Experimental analysis of UPM is both extremely dif- 
ficult and costly to investigate: the need to damp vi- 
brations, assuring adequate vacuum and the necessity 
to use AFM or STM to visualize the results being 
the most important obstacles. 

High costs of experimental setups have led researchers 
to explore the possibilities of computer simulation 
of ultraprecision machining processes. Molecular dy- 
namics (MD) is the traditional method of choice be- 
cause of its relative simplicity and good scalability. 
The increase of available computational power of- 
fers today the opportunity to perform MD simula- 
tions of UPM using 10^-10^ atoms with timescales 
in the order of a few nanoseconds [1], [2], [3]. 

Such simulations are performed in the non-equilib- 
rium (NEMD) regime, with the tool atoms usually 
having their velocity reset to a constant value at each 
step. The workmaterial needs artificial "clamping" 
to ensure it remains in place and does not trans- 
late upon contact with the tool. Periodic boundary 
conditions may [3] or may not [1] be used along 
the direction perpendicular to that of the machining. 
The potentials used are either pairwise, such as Morse 
potential [2] or many-body, such as the Sutton-Chen 
potential [31. Since these are parametrized to repro- 



duce the behaviour of atoms close to their equilib- 
rium positions, they are rather unreliable in situa- 
tions where bonds are ruptured and such is the case 
in the part of the system where the tool tip enters 
the material. In fact it was shown that for brittle frac- 
ture of silicon the widely used Stillinger- Weber MD 
potential gives incorrect predictions, as compared to 
DFT-based methods [6]. 

Apart from UPM there are other interesting sys- 
tems and processes which also do not lend themselves 
easily to experimental investigation and contain re- 
gions where bond-breaking takes place, which renders 
MD simulations thereof unreliable. Among these are: 
metallic nanowires breaking under stress, nanoinden- 
tation, nano-scratching, fracture of nanoscale mate- 
rials. If the system is small (less than 10'^ atoms) 
quantum-based methods, such as DPT or tight- 
binding (TB) which deal with the electronic struc- 
ture explicitly may be used. For larger systems the 
computational cost of these methods quickly becomes 
prohibitive. There is, however, a great need for meth- 
ods that would allow for the treatment of such sys- 
tems, thus allowing to simulate more complex mate- 
rials, auxctics being one of the greatest challenges. 

One solution of this problem is offered by the so-called 
cross-scaling methods, in which a part of the sys- 
tem that needs quantum-based treatment is identi- 
fied and treated accordingly, while the remaining part 
of the system is simulated using the MD method. 
This approach attempts to benefit from the advan- 
tages of both methods, however it suffers from its own 
problems, such as difficulty in determining the posi- 
tion of the region and dealing with boundaries be- 
tween the two partitions of the system. Handshaking 
the two methodologies is the most enduring task. 



The paper is organized as follows. In section 2 we de- 
scribe the formulation of the tight-binding method 
and the molecular dynamics method, later focusing 
on bridging of the two methods. Our computer code 
implementing the cross-scaling is then characterized, 
along with the simulation procedure. In section 3 
we present results of preliminary tests used to vali- 
date the method. Section 4 contains conclusions and 
final remarks. 

2 Computational method 

2.1 The NRL-TB method 

Tight-binding is the only quantum-based method 
that allows for the treatment of systems compris- 
ing several hundred atoms on a typical worksta- 
tion. While relatively simple compared to DFT- 
based methods, it still captures the relevant electronic 
efects, at least in a qualitative manner. The method 
reduces the Schrodinger equation to a generalized 
eigenvalue problem by expanding one-electron wave- 
functions as linear combinations of atomic orbitals. 
Following Slater and Koster [7] matrix elements of the 
hamiltonian are separated into an angular-dependent 
part and distance-dependent two-center integrals, 
which are then parametrized. Brute-force matrix di- 
agonalization which is involved in solving the eigen- 
value problem scales as 0{N^) and is usually the bot- 
tleneck of the calculations. 0(A^)-scaling methods 
have also been developed, but since they rely on bond 
locality, they are less transferable and perform poorly 
for metals [8] . In the tight-binding formalism the en- 
ergy of the system is written as: 

ETB = 2Y,^n + F[n{r)], (1) 

n 

where the first term, involving summation of the elec- 
tronic eigenvalues over all occupied levels 
is the band structure energy. The second term 
is an unknown functional of the charge density n(r) 
and accounts for the remaining density-functional en- 
ergy (coulombic core-core repulsion, canceling out 
the double counting of electron-electron interac- 
tions, etc.) This second term is usually approximated 
by a parametrized repulsive pair potential, which ne- 
cessitates using extra parameters apart from the ones 
describing the two-center integrals. 

In the Naval Research Laboratory tight-binding 
(NRL-TB) total-energy method [9]- [23], developed 
by Mehl and Papaconstantopoulos the second term 
is rid of by clever shifting of the underlying DF 
Kohn-Sham potential. This narrows down the param- 
eter space at the price of introducing a dependence 
of the on-site matrix elements on the local environ- 
ment. It can be shown [9]-[ll] that such approach 
is well-justified and is a generalization of the pair 
potential typically used for F\n{r)\. The NRL-TB 



method offers transferability superior to other TB 
variants and parameter sets for many d-metals are 
readily available. These advantages have made it our 
method of choice for the quantum-based region 
of the system, despite the fact that the computer code 
implementing the method is no longer given away 
by the authors. 

2.2 The molecular dynamics method 

In the MD method a system of classical particles 
interacting with one another via an empirical po- 
tential is followed in discrete timesteps. Analytical 
differentiation of the potential yields forces acting 
on the particles and numerical integration of Newto- 
nian equations of motion gives positions of the parti- 
cles in subsequent timesteps. Aided with the link-cell 
method, it scales as 0{N) and simulations of hun- 
dreds of thousands of particles for timescales of a few 
nanoseconds are feasible on a typical workstation. 
A usual timestep is in the order of 1 fs, we have chosen 
= 2.5 fs. We have decided to choose a fourth-order 
Gear predictor-corrector algorithm as the numerical 
integrator. 

2.3 nanoTB code 

After careful investigation of available non-commer- 
cial MD codes we have decided that none of them 
met our needs for simulating ultraprecision machin- 
ing or other systems involving plastically deforming 
metals. However a parallel code suitable for perform- 
ing NEMD simulations was being developed in our 
group as part of research on carbon nanotubes [4]. 
This code was easily adapted to our task of inter- 
est by the introduction of the Sutton-Chen many- 
body potential. The nanoMD program implements 
the molecular dynamics method and is geared to- 
wards simulations involving external forces, therefore 
making it easy to fix atoms in position, apply con- 
stant translatory or rotational forces, etc. A selection 
of potentials, including Morse, Sutton-Chen, Brenner, 
Finnis-Sinclair, harmonic and anharmonic is avail- 
able. Nose, Nose-Hoover and gaussian thermostats 
have also been implemented. It has been used with 
success for the MD simulation of UPM [3] . 

To perform cross-scaling simulations we have imple- 
mented the NRL-TB method from scratch, incorpo- 
rating it into our code. The program, now called 
nanoTB, accepts the TB parameter files available 
from the NRL website. Two modes of operation are 
possible - in the first a TBMD simulation is per- 
formed, that is the whole system is treated with 
the TB method, and the quantum-based forces it sup- 
plies drive the MD simulation. The second mode al- 
lows for a cross-scaling TB-I-MD simulation - regions 
of the system are designated for which a TB calcula- 
tion takes place, in these regions the quantum-based 
forces replace the forces obtained from MD, while 



the rest of the system is treated with MD. Forces 
are obtained from the TB hamiltonian matrix ele- 
ment derivatives by means of Hellmann-Feynman the- 
orem. Each region may take a spherical, cylindrical 
or cuboid shape and be either fixed in size or dynami- 
cally adjust to accommodate a fixed number of atoms. 
For the cuboid and cylindrical region shapes peri- 
odic or fixed boundary may be used. Charge self- 
consistency may be achieved by means of Hubbard- 
U method or be neglected. Fermi broadening of en- 
ergy levels is used. The regions may be fixed in place, 
move with uniform velocity (eg. along with the tool 
tip) or be made to follow energy peaks in the system 
(useful for tracing defects). 

2.4 Handshaking TB and MD 

Devising a physically sound method for embedding 
the ab-initio region within an MD system is of ut- 
most concern for cross-scaling simulations. Spuri- 
ous effects resulting from the cessation of bonds at 
the TB/MD interface as the TB region is isolated 
from the rest of the system have plagued cross-scaling 
simulations from their inception [5], [26], [27]. The ap- 
proach usually taken to account for the unsaturated 
valencies is to surround the quantum region with vir- 
tual atoms, usually monovalent, which serve to ter- 
minate the dangling bonds [5], [27]. Since they intro- 
duce extra atoms into the calculation, these link-atom 
methods (LAM) increase the computational load and 
require a parametrization of the method not only 
for the species of interest but also for the interactions 
between it and the link-atoms. Even more important 
is the fact that they rely on bond locality which pro- 
hibits their use for metallic systems. In fact it seems 
that no handshaking method that could deal with de- 
localized bonding has ever been proposed. 

Our first attempt to approach this problem was as fol- 
lows. We have conceded to the fact that untcrmi- 
nated valencies will be present in the TB region but 
have instead concentrated on reducing their influence 
on the atomic trajectories. Assuming that the great- 
est disturbances in the TB-generated forces would 
be present in the outer part of the quantum re- 
gion, we decide to give more credence to the MD 
forces for the atoms close to the region perimeter. 
This was achieved by taking for the force acting on 
atom i within the region a weighted average of the TB 
force and the MD force, the weight being the relative 
distance d of the atom from the region centre: 

-d-i^*^^ + (l-d)i^^^. (2) 

This procedure was meant to result in a gradual 
switching from the TB forces to the MD forces 
as the former became more unreliable. 

Unfortunately this simple approach has led to seri- 
ous unphysicalities in the test cases we have studied. 
Upon centering a spherical TB region of 14 A in diam- 



eter on a self-interstitial defect in Cu, the quantum- 
based forces acting on the atoms close to the perime- 
ter of the region were two to three orders of mag- 
nitude greater than the corresponding MD forces. 
This was due to the fact that these atoms were 
only very slightly disturbed from their lattice posi- 
tions, yet they were heavily influenced by the isola- 
tion of the TB region from the rest of the system. 
As a consequence, even with d = 0.95 the force com- 
puted from (2) was overestimated by an order of mag- 
nitude or so for all atoms close to the region bound- 
ary. A similar scenario occured for a cylindrical region 
centered on one of the edges of a cuboid of 1638 cop- 
per atoms, mimicking a workmaterial to be machined. 
The strong TB forces resulting from the isolation 
of the region easily outweighed the relatively small 
MD forces acting on these same atoms. As the hamil- 
tonian resulting from (2) is not conservative, this re- 
sulted in constant increase of system energy, effec- 
tively melting the surrounding material. 

Thus it became obvious that a more sophisti- 
cated method of combining the TB and MD forces 
was needed. We have settled for a method similar 
in concept to the one described, but using a more 
aggressive mixing in of the MD forces. We have de- 
cided to cut-off the influence of the TB forces using 
a smooth function similar to the cutoff function used 
in the NRL-TB method [12]. Instead of (2) we use 

Fi=w{d)F,'^'' +w{l-d)Fr, (3) 

where 

w(d) = 5 -, (4) 

^ ' l + exp(-f +5)' ^ ' 

with I = 0.06. This leads to a non-linear weighing 
of the TB and MD forces. 
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Fig. 1: The shape of the weighing function w{d). 
Note u;(0) « 0, w{l) w 1. 

This improved embedding has proven more success- 
ful against our test cases. As an example consider 
a cylindrical quantum region of diameter 13.25 A 
placed on a surface of a cuboid, with the cylinder 




axis parallel to the cuboid edge. This system consist- 
ing of 4305 total atoms and 120 TB atoms was simu- 
lated for 8 ps and results of the two mixing methods 
were compared. 
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Fig. 2: Artifacts introduced by unterminated valen- 
cies of the TB region using the (a) linear; (b) nonlin- 
ear mixing method. A snapshot of the system config- 
uration at t = 0.85 ps is shown, with the TB atoms 
drawn in darker colour. 



The artifacts introduced by the nonlinear method 
were markedly smaller, although still visibly present. 
This, naturally, comes at the cost of having to ex- 
tend the region in practical applications as the atoms 
further than d = 0.5 from the region centre are effec- 
tively driven by the MD method, however their inclu- 
sion in the TB calculation serves to create a proper 
local quantum environment for the atoms closer to 
the region centre. 

2.5 Simulation procedure 

As a preliminary test of our code and the non- 
linear mixing method, we have performed a cross- 
scaled simulation of nanoscratching of a copper 
workmaterial with a copper tool. The workmaterial 
was a cuboid of 20x10x5 fee unit cells, for a total 
of 4305 atoms. The tool consisted of a cuboid of 8x8x5 
fee unit cells (1445 atoms), rotated by 45 degrees 
along the z axis and moving along the [010] direc- 
tion. Periodic boundary conditions were used along 
the z (in-plane) direction. Atoms on the bottom sur- 
face of the workmaterial were fixed to prevent it from 
translating upon contact with the tool. For simplicity 
the tool was assumed to be infinitely hard (all forces 
acting upon it were ignored.) the tool moved with 
a uniform velocity of 10 m/s, which is considerably 
more realistic than in other simulations of UPM [1] , 
[2]. A pure MD simulation using the Sutton- Chen po- 
tential was performed for comparison. The quantum 
region was cylindrical with periodic boundary condi- 



tions along the z axis and was centered on the tool 
tip, moving with the same velocity as the tool. It was 
taken to be only 8 A in diameter for efficiency rea- 
sons, in realistic simulations the region should have 
a radius of at least 6.62 A, since this is the cut- 
off distance for Cu in the NRL-TB method. Since 
the region was fixed in size, the number of atoms 
contained within it varied with time, with an aver- 
age of 70 atoms. The simulations were performed at 
300 K with a Nose-Hoover thermostat having a time 
constant of 250 fs. A timestep of 2.5 fs was used. Ini- 
tially the tool tip was at a distance of 4.75 A from 
the workmaterial, giving it ample time to equilibrate 
before contact was made. 



3 Results and discussion 

We consider the presented simulation to be merely 
a preliminary test, therefore we look at the obtained 
results only qualitatively. The differences between 
the cross-scaled and the reference, pure MD simula- 
tion are rather subtle, but noticeable. For both sim- 
ulations we observe the jump-to-contact (JC) phe- 
nomenon when the distance between the tool tip and 
the workmaterial is about 4 A. This is in agreement 
with pure MD simulations described in [1], [24], [25]. 
Since the tool is infinitely rigid, it is the workmaterial 
that deforms slightly to contact the tool. At this point 
the cross-scaling and the reference simulation are ex- 
pected to perform identically, because as the TB re- 
gion is centered on the tool tip, at this distance there 
are only 15 workmaterial atoms within the region and 
these are driven almost exclusively by MD. 
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Fig. 3: The jump-to-contact phenomenon. The TB re- 
gion is marked with a circle with the TB atoms drawn 
in a darker colour. Note how this effectively reduces 
to an MD simulation as the workmaterial atoms have 
only begun to enter the quantum region. 
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Fig. 4: A snapshot of the configuration at t = 100 ps, 
after the tool has penetrated the workmaterial 5.25 A 
deep, a) cross-scaling TB+MD simulation (TB atoms 
shown in darker colour); b) pure MD simulation. 



Until the distance between the tool tip and the work- 
material reaches about 1.5 A, the values for the nor- 
mal force acting on the tool are practically identical 
in the two simulations and extremely small - the tool 
is attracted by the workmaterial with a force lower 
than 0.1 nN. This force fluctuates quasi-periodically 
as the workmaterial "wiggles" under the tool in re- 
sponse to the attraction. As the tool tip gets closer 
than 1.5 A to the surface, the force becomes repelling 
and then rises steadily after contact as the tool pen- 
etrates the workmaterial. The forces achieve their 
highest value of about 0.8 nN at the penetration 
depth of 8 A (where the simulation terminates) . Dur- 
ing the course of the simulation, as more atoms are 
treated with TB, the differences between the simula- 
tions become noticeable. The normal forces observed 
in the cross-scaling simulation are seen to be 10-30% 



greater than the corresponding pure MD forces. 

Taking a closer look at the configuration of the sys- 
tem after the tool has entered the workmaterial 
5.25 A deep, we note slight differences between 
the cross-scaling and reference simulation. A compar- 
ison of the figures below reveals that the cross-scaled 
simulation shows a smaller disturbance of workma- 
terial several layers under the tool tip. Slight differ- 
ences in the pressure maps (Fig. 5) and the fact that 
the normal force experienced by the tool is greater 
for the cross-scaling simulation lead to a supposition 
that the material is in fact slightly harder than pure 
MD simulations would predict. 
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Fig. 5: Pressure map (arbitrary units) for the snap- 
shot in Fig. 4. a) cross-scaling TB-I-MD simulation; 
b) pure MD simulation; c) the difference between 
the two (image of system configuration overlapped 
for clarity). 



4 Conclusion 

We have developed a computer code that allows for 
cross-scaling simulations using the methods of NRL 
total energy tight-binding and molecular dynamics 



in combination. The program allows for the embed- 
ding in an NEMD simulation a number of quantum- 
based regions that are dynamic in size and position 
and may take various shapes. We have considered two 
methods that would tentatively deal with the prob- 
lem of unterminated valencies in metallic systems 
on the TB /MD boundary, showing the first one to fail 
and the second to be moderately successful. The non- 
linear mixing method minimizes the surface effects 
on region boundaries at the price of having to extend 
the region, which increases the computational load. 
Also the hamiltonian resulting from the application 
of the method is not conservative. We have tested 
the code and the nonlinear cross-scaling method 
on several simple systems, finding the code to work 
and the method to behave better than the linear 
one we have proposed earlier. For the simulation 
of nanoscratching, our cross-scaling method gave sim- 
ilar, but not identical results to the pure MD simula- 
tion, we bear in mind, however, that the region size 
used for the test was too small for practical applica- 
tions. We envisage adapting the recent learn-on-the- 
fly method [28] to the Sutton-Chen potential in order 
to attack the problem of embedding quantum-based 
calculations within an MD framework from a different 
perspective. 
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